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ABSTRACT 


Linear prediction has become an important tool for stationary time series analysis. 
The all-pole system model which is a product of the linear predictive approach has ap- 
plications in numerous engineering problems. This thesis develops a simple method for 
obtaining the two-dimensional all-pole system model. The lattice structures that can be 
used to implement the prediction error and synthesis filters are also shown to have an 
analogous two-dimensional counterpart. The construction of these filters is in terms of 
orthogonal Szego polynomials which can be used to solve the two-dimensional block 
Toeplitz normal equations in a recursive manner. This recursion not onlv leads to the 
two-dimensional (2-D) lattice structure, but also allows for expansion of the filter order 
without resolving the normal equations. Several examples are presented using the two- 


dimensional linear prediction results for spectral estimation and signal synthesis. 
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I. INTRODUCTION 


A. LINEAR PREDICTION 

The initial goal of linear prediction is to estimate the present value of a random se- 
quence based on a linear combination of all the past values. This is equivalent to de- 
termining a linear, shift invariant, causal prediction error filter which whitens the 
random sequence. There has been a great deal of study of one-dimensional linear pre- 
diction analysis and the usefulness of the resulting filters. The purpose of this thesis 1s 
to extend in a straightforward manner the results of one-dimensional linear prediction 
to two-dimensional signals. 

There are several efficient methods for the solution of the normal equations that 
arise in one-dimensional linear prediction problems. The most familiar of these methods 
is the Levinson algorithm. The Levinson algorithm alleviates the need to perform a 
matrix inversion in the solution of the normal equations for the filter coefficients. The 
recursive nature of the algorithm also allows the prediction error filter to be implemented 
in a lattice structure. These two features will be of primary concern in the solution of 
the two dimensional problem. 

There are two general approaches to linear prediction. The first method, which 1s 
usually referred to as the spectral factorization method, begins with the knowledge of the 
power spectral density of the random process and attempts to find the model parameters 
that fit the spectrum. This works well in one-dimensional applications but cannot be 
extended to two-dimensional signals. In general, a two-dimensional spectrum is not 
factorable into two polynomials. Although progress has been made in factoring two- 
dimensional spectra in terms of the complex cepstrum and the results are theoretically 
exact, these theoretical results cannot be attained in practice. We will, therefore, focus 
our attention on a second method of attaining the prediction error filter. 

The solution method dealt with in this thesis begins with the individual samples of 
the random process and develops the prediction-error filter from these samples. This 
method is Known as autoregressive model fitting or all-pole modeling. The solution will 
be developed in terms of orthogonal polynomials that will not only yield a Levinson-like 


algorithm but are easily extended to two-dimensional signals. 


B. ONE-DIMENSIONAL LINEAR PREDICTION 
1. All-Pole Model 
We are given a finite length sequence of the past values of the random sequence 
Y(n). The problem is to estimate Y(n) based on a linear combination of these samples. 
For the simplest case We can estimate Y(n) based solely on the previous sample. The 


estimate, denoted by Y(n), is 
¥(n) = aY¥(n—1) (1) 


The prediction will have some error, e(”), which will be minimized in some sense 
by the choice of a. The least-squares minimization criterion is commonly used. This 
criterion provides a simple approach to the solution of the linear-prediction problem and 
thus the modeling problem. We proceed as follows. 


First. we redefine the estimate as 


N 


Yin) = -aY(n—-—1) . (2) 
Now the prediction eriorm yaliebe 


¥(n) — Y{(n) 
Y(n) + a¥(n—-1) . 


e(n) 3) 


| 


The squared error is to be nunimum with respect to the coefficient a. This is 
easily accomplished using derivatives as shown below. Differentiating [e?()] with re- 


Spect to a gives 


£ [e*(u)] = 2e, He, = 0 (4 
or 
r¥(n) + aY(n — WILY) + aVqeieMe (5) 
which yields 
[¥Ga) sar Ga 7 — le (6) 


Finally, solving for a 


Y(n) Y(n — 1) 
SSeS ) 


Since the value of a must minimize the error, e(7), for all values of n, the value 


of a in (7) could be rewritten as 


n 


> YOYG=1) 


oe ees en. (8) 


(=a 


i=] 


Note that (8) will select a such that }°e?(7) is minimum. It can be seen that the 
coefficient a is a function of the autocorrelation elements of the given sample space, and 
this method 1s usually referred to as the autocorrelation or Yule-Walker approach. 


If the order of the predictor is now increased such that 


+ 


= —[a¥n—-1l) + a ¥(n-2) +. | (9) 


the prediction error becomes 


A 


wi Yin) = Veni aii — 1) 2 a — 2). (10) 


If (10) is considered to be the difference equation of an all-zero filter, the 
prediction-error filter can be described by taking the z transform of (10) as 
eee 2 )Alz) (11) 
where 


A(z) = 1+ a,z7' =H ay2~ ee (12) 


Now if (10) is rewritten and e() is taken to be the input to the system, the es- 


timated signal generator model results 





where 


I 
B(z) = = er (14) 
Do 2 ez ees 


It can be shown that e, 1s a White noise sequence with variance E[e?]. This im- 


plies that Y(z) and, consequently, Y(n) can be generated from a white noise sequence 


as shown in Figure |. 





Figure 1. Prediction error and system model filters. 


This assumes that the sequence A(z) is causal and causally invertible. If this 1s the case, 


we can generate the error sequence from Y(n) and vice-versa by an invertible filter op- 


eration [Ref. I]. The signal generator model, B(z), is the best estimator in a least-squares 
sense and A(z) 1s the prediction-error filter. These filters can be used in a number of ways 
and the simple solution for the filter coefficients is the primary reason for the popularity 
of linear prediction in signal processing. 
2. Normal Equations and Orthogonality 

Recalling that we wish to minimize the sum of squares of the error sequence 
leads directly to an important result in linear prediction known as the orthogonality re- 
lationship between the sample space and the error. This principle states that the error 
e(m) 1s orthogonal to all the samples of Y(n) used in obtaining the estimate Y(n). That 


1S, 
Peeler) (a — a1) nl — Oey = oy) 


where N is the order of the predictor. This result 1s derived for the two-dimensional case 
in Chapter If. It is the central idea behind the solution of the linear prediction problem 
in terms of orthogonal polynomials. It also leads directly to the normal equations in 
both one- and two-dimensional formulations. 

The derivation of the normal equations is a simple extension of the 
orthogonality relationship (17). Suppose we seek a kth-order predictor and the value 
of the prediction-error variance. We will require (k+1) linearly independent equations 
to solve for the k prediction coefficients and the prediction-error variance. These coef- 


ficients can be derived as follows. 


Y(n) - (71) 
Yn) + a,¥n—1) + @¥(n—2) +--+) +. a,¥(n—k) 


Cn 
(15) 


Substituting (18) into (17) gives 
EL¥(a) + a ¥(n—-1) + @¥(n-2) + = ,¥(n-)D] =O LSI N (19) 


If we now define the coefficient a,=1 and note that for a stationary process the 


mrocoireldation, Aix), 1s given bY EL Y(R)Y() | = "R,_,, (19) may now be written as 


N 


El) a¥(n—k). Y(n—i) | = 0 (20) 


k=O 


N 
> GRi-k) = 0 sis (21) 
k=0 


where N 1s the order of the predictor being used and R is the autocorrelation function 
of the random sequence Y(n). The N equations gencrated by (21) can be written in 


matrix form for N=2 as 
Ay 
ay == 5 (22) 
ee IR, 0 


The error variance € 1s given by 


c= Le] = Elen Ya) — Yo] 


ae (23) 
=: Eleeaie ale 
Substituting (18) into (23) yields 
N 
c= E] Ya ¥n— 4), Yn) (24) 
k=0 


P 

Which in terms of the autocorrelation becomes e = > a,R_,. For example, with N= 2, 
k=0 

this can be combined with (22) to vield 


Ry Ro areal 
Ry oi meee 
R, R, Ko Qy 


(25) 


fy 
i 
Ze, 


This autocorrelation matrix is seen to be (N+1) by (N+1), Toeplitz, and since 
R_, = K,, itis svinmetric. The various solutions for (25) occupy much of tne literature 
on linear prediction. For our purposes, the desired solution must provide a simple cx- 
tension to tWo-dimensional signals and be adaptable to a lattice structure. J.H. Justice 
(Ref. 2} has derived such a solution in terms of Szego polynomials. Although quite 
similar to the Levinson algorithm, it has the unportant advantage of not relving on the 
symmetry of the correlation matrix. The two-dimensional autocorrclation matnx will 


be found to be Toeplitz in blocks, and, although this provides some simplification, the 


symmetry of (25) is lost. Fortunately, a solution in terms of orthogonal Szego 
polynomials is easily derived [Ref. 3]. 
3. Applications of the Autoregressive Model 

Linear prediction plays an integral role in the solution of a number of engi- 
neering problems. The resulting model can be used for high resolution spectral esti- 
mation [Ref. 4]. Speech signals can be reduced from a large volume of data to a 
comparatively small data space by linear predictive coding [Ref. 5]. A model of the 
generating function for the random sequence can be derived. Lattice filters can be im- 
plemented which can be used to form the prediction-error filter in the forward direction 
or synthesize the random process in the inverse direction. 

The many applications of linear prediction in time series analysis has generated 
a great deal of interest in extending the results to multi-dimensional random processes. 
Multi-dimensional signals occur in situations where the samples are spatially dependent 
such as radar, sonar, and image processing. The signals considered in this thesis are wide 
sense stationary and scalar-valued sequences of two variables. The two-dimensional 
autoregressive model that is derived will be used in some of the same applications clas- 
sically associated with the 1-D model including spectral estimation and sequence gener- 
ator estimation. An alternative implementation of the model in a lattice structure will 


also be presented. 


C. NOTATIONAL CONVENTIONS 
I. Mfatrix Notation 
This thesis will deal with a finite size rectangular sample space. The individual 
samples will be denoted by ¥(m, —A,n,—-, where & and / are measures of distance from 
the point at which the value of the sequence is being estimated. The sample space will 
contain Nl] x N2-1 samples. Figure 2 shows an example of the sample space with 
N1=N2=4. In this case we would be predicting the value of Y(3,3) based on the 15 
samples available. Notice that, in this case, A and / have a maximum value of 3. Spa- 
tial samples are indexed by nl and n2 to denote rows and columns, or the vertical and 
horizontal directions, respectively. 
2. Causal Samples 
In one-dimensional signals the clear meaning of past and future leads to a 
Straightforward definition of causalitv. In two-dimensional signals the ordering of the 
points is of central importance; therefore a reference for past and future needs to be es- 


tablished. In this thesis causal samples will be taken to be in the quarter plane shown 





Figure 2. Sample space for NI =N2=4. 


in Figure 3. This 1s equivalent to describing the desired filter as a third-quadrant filter. 
Due to the symmetry of the 2-D autocorrelation the third and first quadrants produce 
the same results. Hlowever, causal will be used only m reference to the samples shown 
1 Figure 3. 
3. Polynomial Notation 

Two-dimensional polynomials will be denoted by P,, where & 1s the degree of the 
polynonual m the nl direction and /1s the degree in the n2 direction. When the need to 
refer to an individual element in a polynomial arises, the element will be tdentified as 
Pity), Where ¢ and s give the power of the bivariate variable. The correspondence be- 
tween the z transform of the system model and the coefficients in the orthogonal 
polynomials leads to the use of z, and z, as the independent variables tn this study. 


Thus the polynomial P,, will be a 3x3 bivariate polynomial consisting of the elements 


a 


____ Causal 
eee id 





Figure 3. Causal Data Space 


P27(0,0) P22(9,1)z, 
P= P2102, po2(1, 12,2, 
2 2 
P22(2,0)2) P22(2,1)27 22 


Y(n1,r2) 


; (26) 
2: 
A 


The development of P,, is in terms of positive powers of z. The desired filter 


must be causal and therefore in negative powers of z. This is easily accomplished by 


matiplving P,, by zj*z;'. 


The resulting causal filter will be denoted A,, . 


The leading 


coeflicient of P,, is p,({k,J) which will become a,{0,0) in the causal filter. 


Il. TWO-DIMENSIONAL LINEAR PREDICTION 


A. PROBLEM DEFINITION 

The two-dimensional linear prediction problem 1s conceptually identical to the one- 
dimensional problem. That is, we attempt to predict the value of Y(nl,n2) based on a 
linear combination of the previous NIXN2 rectangular block of data. The coefficients 
which result in the optimum prediction will comprise a linear, shift invariant, causal, 
prediction error filter. The problem is to find the coefficients that minimize in some 
sense the error between the actual value of Y(n1,n2) and the predicted value, ¥(n1,n2). 
This will require the solution of a block-Toeplitz matrix which contains NIXN1 blocks 
and each block will contain N2xN2 correlation coefficients. We will require the solution 
to provide the means for a direct realization of the filter or implementation in a lattice 
structure. The mean square error criterion 1s again utilized to obtain the optimum pre- 
diction filter. If a solution to the desired filter can be achieved recursively, then the re- 
lationship between successively higher order filters can be used to implement the lattice 
filter. 

We begin bv seeking a solution to the linear prediction problem, and will assume 
that Y(nI.n2) is a two-dimensional stationary random sequence. Given the causal block 
of NIXN2 data samples we need to generate ¥(n1,n2) from a linear combination of the 


data that minimizes the prediction error. The prediction, ¥(nl,n2), can be written as 


A 


V(al2) = — ayy V(r — 1,22) —a9,(m1,m2 — 1) — a, Yl — 1,v2—-1)— -- 07) 
— Oni ho) ¥(n] = AI + ley a NZ + 1) 


where the a,,'s are the filter coefficients associated with data samples that are delayed by 
sin the nI direction and r in the n2 direction. The coefficients can be assembled in a 
matrix to yield the desired two-dimensional linear prediction filter. This filter will be an 
all-zero model with the finite sample space Y as the input and the prediction, Y, as the 


output as described below. 


1Q 


oNT=,N2—) °° Ani~1,) 4N1~1,0 


y= Y 2-D filter coefficients (28) 


Q) N2-1 Ae ay) aio 


Qo N2—1 ee Qo) 0 





Figure 4. Generation of the prediction error 


The prediction process can be depicted as in Figure 4. The polynomial A(z, , z,) 1s 


simply the z transform of (27) given bv 


~1 ~1 —1_-2 
Al2y 522) = —(4i92) + A122 $42) 220 + oe) (29) 


For any choice of the coefficients some error will be present. This prediction error is as 


seen in Figure 4 and can be written as 


Ni-1 N2—1 
e=¥(nl.n2) + YY ay ¥(nl —syn2—1) « (5, # (0,0) (30) 


5—0 9 =U 
We will define the filter coefficient a, = 1 and condense (30) to the form 


N1—-1 N2-1 
e= > >) ay¥(nl—syn2—2) | (31) 


s=0 r=0 
The task now 1s to find the filter coefficients that minimize the error. 


B. TWO-DIMENSIONAL ORTHOGONALITY 
1. Description 
The error, Y(nl,n2)- Y(nl,n2), can be visualized as in Figure 5. Any linear 
combination of Y(nl-1,n2) and Y(n1,n2-1) will vield an estimate of Y which lies in the 
plane of Y(nl-1.n2) and Y(nl.n2-1). If Y(n1,n2) lies in this plane, then Y will be exact 
and there will be no error. For any other Y(n1l,n2) the error will be minimum if it 1s 
orthogonal to the sample space spanned by Y(n1-1,n2) and Y(n1,n2-1). This is equiy- 
alent to defining ¥ as that portion of Y(n1,n2) that is parallel to the sample space. Thus 
Y-Y leaves onlv the component of Y that is perpendicular to the sample space as the 
error. Although it can not be drawn for higher order predictors, the minimum error 1s 
always such that it is orthogonal to the space spanned by all the samples of Y used in 
the estimate. 
2. Derivation 
The orthogonality of the error to the sample space can be easily derived. First 


we define the inner product, <X,Y>, of two sequences, X and Y, by 


<4 = >) ae (32) 
k l 


Now the inner product of the squared error can be expressed as 
2 172 
<e°> = <(Y—Y)*> (33) 


and we seck a minimum squared error with respect to the filter coefficients a. 
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Figure 5. A visualization of tle orthogonality principle for two samples. 


ee) (34) 
da a 
Or 
pac nas) (35) 
da 


Where 7 are the partial derivatives with respect to the elements in the coefficient ma- 
are; 


trix. Substituting (30) for e 





Mie Ne| 
d dares : 
meme aa | ee et a, ¥(nl — s,n2—1 36 
7 ‘ Dal ) 3 » st¥( J (36) 


leads to 


Se = (Hn — 1,72), ¥(n1,n2 — 1), , Yl — N1 = 1,02 — N2 = 1)]) ee 
dQ 


Substituting (37) into (35) yields 


<2, Y nl, — sl — ty — Cn se al 
Oe N2 (38) 
hy on 


This result, (38), states that the inner product of the error and each of the past 
samples of Y(nl,n2) used in the estimate is zero which requires that the error be 
orthogonal to the sample space. Not only is the error orthogonal to the sample space, 
but errors created by successively higher-order predictors are orthogonal to each other. 
Also, the sequence e(n,, 7,) will be uncorrelated to the degree that the order of the pre- 
dictor allows. That is, if we had access to an infinite causal data space, then we would 
produce a purely white error sequence. Since we will always be limited to some finite 
size filter, we will produce an error sequence that is partially correlated. 

3. Two-Dimensional Normal Equations 

The sequences considered in this thesis are taken to be wide-sense stationary. 
If thev are also taken to be zero-mean signals, the inner product of two sequences 1s 
equivalent to the expected value of the product of the sequences. The expected value 
of Y(nl,n2) and Y(nI-s,n2-t) reduces to a two-dimensional autocorrelation that is solely 


a function of the delavs s and ¢. Thus, 
< Y(nal,n2),¥(nl —sn2—y)> = ELY(nl,n2),¥(nl —s,n2—1)] = R,, (39) 


where R will be used to denote the autocorrelation with lags s in the nl direction and 
in the n2 direction. Equation (38) can now be used to form a block-Toeplitz matrix that 
is analogous to the Toeplitz matrix generated by the normal equations in one- 


dimensional linear prediction. Thus, we can write 


Ele, ¥Ynl —s,n2—t)]=0  (s,t) # (0,0) (40) 
substituting for e 
MEIN 
E} Y(nln2)+ > ay¥(nl — kyn2—D, Mal — sn2—-) | = 0 O<s<N1 
k=0 10 0<1<N2 Tam 


(s,£) + (0,0) 


that can be written as 


mils WS) O<s< Nl 
a,R(s—kt—-)D =0 .0st<N2 (42) 
k=0 /=0 CHE 


For example, with a N1l = N2=2, and with the first filter coefficient defined as 


Qj, = 1 the normal equations become 


E\ ¥(ni,n2) + ay, Y(n1,n2 — 1) + ayy Y(nl — 1,2) + 


(43) 
+a), ¥(nl —1,72—1), Yal—sn2—)| = 0 


where (42) or (43) are the two dimensional normal equations. They are a direct result 
of the orthogonality relationship described by (40). These equations can be written out 


to yield NIxN2-1 linearly independent equations. For example, with NI=N2=2, we 


have 
S=O7=1 Ro, +49) Rog + Ay9R_1 1 + 4, R_19 = O (44) 
S=1,t=O0 Rig + aR, 1 + Ayo Rop + 4), Ap _) = O (45) 
s=lr=1 Ri1 + MH Rig + MoRo) + 4; oo = O (46) 


Equations (44)-(46) can be augmented by the expression for the squared error 


€ Which is given by 


ee ewes 0 (47) 
OT 
Nl=—] N2—1 
= E) > Yay ¥(al — kyn2 — 9, ¥(n1,n2) (48) 
k=0 /=0 


and in terms of the autocorellation 


Pmt 


ein. 


to 


=) 


ayjR(—k,— 2) . (49) 


f 


k=0 


~ 
ll 
=) 


Assembling the above equations into a matrix yields 


Roo Ro -| R_, 0 R_, =) Aog 
Roy Kop Riya R19 4) 
Rio Ry Roo Ro -] Ai9 
Ri, Rig Roy Roo Q)y 


(50) 


ce aes PG 9 


The resulting matrix is (NIN2)x(NIN2) square and will consist of blocks which 
are Toeplitz. There will be NIxN1 blocks and each block will contain N2xN2 corre- 
lation coefficients. The matrix can be seen to be Toeplitz with respect to the blocks and 


ella be 
DO, DM) |] 4, 0 


Where the @,’s are blocks in the autocorrelation matrix and H indicates hermitian 


can be written as 


transposition. Each @, is an N2xN2 Toeplitz matrix. For the example from (50), we 


have 
R R 
i 00 - (52) 
Ko; — Roo 
R R 
Ry; Ryo 
—— (Ago, aaa (54) 
a = (40, 4) (55) 
@ = (e, 0)" (56) 


If an element in D, is R,, , the corresponding element in ©” is K_.. Simgeue 
are dealing only with real data the conjugate notation can be dropped from R. It 1s also 
worthwhile to note that RK, = R_,_,. 

The solution of this matrix for the coefficients a,, will yield the two-dimensional 
whitening filter. There has been a great deal of research involving block Toeplitz ma- 
trices and there are several algorithms that could be used to obtain the coefficients. In 
order to implement the filter in a lattice structure, a recursive solution beginning with a 


zeroth-order filter and proceeding to the desired order is necessary. In one-dimensional 


prediction error filtering, the Levinson algorithm successively represents the higher order 
filters in terms of reflection coefficients. It will be seen that an analogous solution to the 


two-dimensional prediction error filter can be developed in terms of Szego polynomials. 


C. SZEGO POLYNOMIALS 

The solution to the system of equations Ma =e does not require that the inverse of 
the matrix ® be found. All that is required is an orthogonal set of polynomials that span 
the space defined by ®. This result follows the orthogonality relationship of the pre- 
diction error and the sample space. For example, if we had a second order and a first 


order prediction error filter the corresponding prediction errors would be 


e, = YI) —a,(1)¥(0) (57) 
ey = ¥(2)—a,(1) ¥(1) — a2) ¥0) (58) 
Eley, e)] = Ele, ¥(1) ~ a,(1) ¥(0)] (59) 


But (44) requires that e, is orthogonal to both Y(1) and Y(0), and therefore 


Ele,e,J = 0 . (60) 


It is evident from the above observations that the choice of A,., must be such that 
it is orthogonal to A,. If we begin by defining A, as any convenient value and proceed 
to find the successively higher-order orthogonal polynomials until the desired predictor 
order is reached, the final polynomial will hold the coefficients which will provide the 
optimum prediction-error filter. To accomplish this we will use Szego polvnomuals. 
Szego polvnomials are a special class of polvnomials which have zero inner products 
@recran interval on the unit circle [Ref. 6]. This is precisely the requirement of the 
orthogonal polvnomials relating successively higher order-prediction errors. J.H. Justice 


[Ref. 3] has presented a recursive method of generating Szego polvnomials of the form 
Pe = Pat + Ppt’) + + Po (61) 
where 
z= eK (62) 
The matrix @ is an inner product matrix which may be utilized to orthonormalize a 


sequence of polynomials on the unit circle [Ref. 6]. The resulting orthonormal sequence 


7 


wil be composed of Szego polynomials. If the filter sequence A, is determined in terms 
of Szego polynomials that are orthogonal over the space associated with the 
autocorrelation matrix , the prediction error filter will result [Ref. 2]. Further, the Szego 
polynomials can be determined in a recursive manner which is necessary for the repre- 
sentation of the filter in a lattice structure. In a one-dimensional filter this will result in 
the Levinson algorithm. J.H. Justice [Ref. 2] has developed the bivariate extension of 
Szego polynomials that can be utilized to orthonormalize the block Toeplitz matrix 
generated by the two-dimensional autocorrelation of the sequence Y(n1,n2). 
1. Levinson Algorithm in Terms of Szego Polynomials 

Any positive definite matrix can be defined as an inner product space and can 
thus be used to define orthogonal polynomials spanning that space. The autocorrelation 
matrix produces such a space. If the autocorrelation matrix is generated from a Se- 
quence, Y= Y,, Y,,----,¥y, then the sequence of polvmomuals, 1) 2,z*)-9y2 ecu 
orthonormalized with respect to the autocorrelation matrix. The resulting polynomials 
are P(z). P,(z). -, P:{z). The subscript on P denotes the degree of the pols nemuaiauuaac 
requirement on P, is that it is orthogonal to all polynomials of lesser degree. Thus, we 


have 
<P Pie = olh—7) (63) 


A one-dimensional sequence of orthogonal polynomials will consist of the terms 


Po = Poo 
Pi =Piy\Z + Pio 
(64) 
N N=] 
Py = PynZ + Prn—Z ota eee NY 
and in general, 

k 
j=0 





Justice’s method for determining these polynomuals recursively results in an al- 
gorithm equivalent to the Levinson algorithm. This algorithm for the 1-D case is briefly 


summarized as follows. We begin by setting 


Py = Rp? . (66) 


Then the subsequent polynomials can be determined by the recursive relationship 


papauslie ‘ 
Pr4i(Z) a (1 a Ide | *Pan) - CZ@ - Benne nZ)) (67) 
where 
Pr(z) = 2"P,(z~’) (68) 
and 
Sia > ResiPrk (69) 
k=0 


The equation (67) can be expressed 1n the form of a causal prediction-error filter 
that vields the familiar form of the Levinson algorithm. The nth-order polynomial 


produced by (67) will be of the form: 
Pi(2) = Prat” + PaniZ to + Pro (70) 
The causal prediction-error filter is of the form: 
An(2) = yg + Gq yz) be + Aye” (71) 
where we let a, = p,,, and A, = z"P,. Then multiplying (67) by z~” yields 
2 Pai) = 2" Py(2) 2A Pal) (72) 


where the normalizing factor, (1 — |/,|2p2,)-7, has been absorbed into the definition 
@uepeand the leading coefficients are defined as p,, = a, = 1. Now, we define the 


backwards polynomials as follows 


19 


Aj(2) = 2"A,(2~ 
(z) a, (2 ta 74) 
=z (z ea Gs )) 
P,(z) = 2”A,(z) (75) 
where (72) can now be written as 
Anai(2) = A,(z) — 2 '4,A4;(2) (76) 
Angi(2) = 2 7 P(4,(27') — 24,4,(2)) (77) 
= 2 'Ay(z) —2-"4,(2"4,(z)) 
Aj(2z) = 2 An(2) — ApAn(2) (78) 
ae 
}, = = ———— (79) 
ee, 
—— 


Equations (76), (78) and (79) represent the usual Levinson recursion of one- 
dimensional linear prediction problems. The coefficients found using the Levinson al- 
gorithm will be developed in negative powers of z. The coefficients found using (67) are 
in positive z and differ by a normalizing factor. Both sets of polynomials, A and P, will 
be orthogonal sets of polynonuals on the unit circle spanning the same inner product 
space. Since the definition of the 4, is simpler using (67), we will proceed to the two 
dimensional extension using positive powers of z and make the necessary conversion to 
a causal filter after the coefficients have been found. 

2.  2-D Extension 

Now we want to extend these ideas to two-dimensional filters. The necessary 
extension to bivariate orthogonal polynomials has been developed by Justice [Ref. 2]. 
The derived polynomials will be orthogonal in two variables with respect to the inner 


product space of M as described by 


ECP.» Py) = 6(k—sS\o(l—2) . (SO) 


The set of polynomials to be orthonormalized are shown tn Figure 6. The procedure 
Starts with the first row which is a function of z, alone. This row, therefore, can be 


orthonormalized by the same method presented for one-dimensional systems. 


1 2 
1 ZZ y a. . 
2 2 
2 1 1 12 
1 LL, 4a ae 
7 Fy’ 29 
, ia LZ ee e@ ee 


Figure 6. Variables used to form the bivariate orthogonal polynomials 


This 1s equivalent to taking the block @, from the two-dimensional autocorrelation ma- 


tmx and finding the filter coefficients a), doa, + 5 Gam Thus, 


Rog Rony oe Rom 
Roy Rog Romi) 
Dy = en Ail (81) 
R,A0 tes Roo 
where 
_ ip 
ly = (9, gis + Fam) + (82) 


It 1s easily shown that Rk, = R,_,, so that this is a symmetric matrix and a re- 
cursive solution has been presented. Once the first row of orthogonal polynomials has 
been determined, the polynomials of the next row can be expressed as a linear combi- 
nation of the lower-order polynomials. For exainple, with N2 = 3 the first row will have 


polynomials Py), Po, and P,,. The first polynomial in the second row to be found ts 


P\. This can be accomplished as follows: 


21 


First note that if Po is multiphed by z, its order is changed and it is no longer a 


member of the orthogonal set. We may then write 
21Poo = 4ooPio + 4o1P01 + 402/02 (83) 
and after multiplying on both sides by P,, and taking the expected value we find 
Elz; Poo, Piod = AooELP ios Pio] + 401£LPo1, Piol + Zo2EL Poa, Piol - (84) 


Since P,, 1s orthogonal to all polynomials of lower order in z,, this reduces to 


Elz; Poo: Pio] = Aoo (85) 
and similarly 
Elz; Poo, Poil = 4o1 (86) 
Elz; Poo. Pool = 4o2 - (87) 
Thus, 1n general we find that 
k N2=1 
UE DePha YiaPae in d Da (88) 
s=0 [= 


Now we define the reverse polynomial 


P21, 22) = ade 9 <9 Se Ora =k 
O<j<em (89) 


m= N2—] 


Note that P: will also be an orthogonal set of polynomials and will span the same space 


as P,. Now (88) may be written in terms of the reverse polynomials 


Z) Pp, = ~~ At ear ar Sat k+1,1 ar S Sab Ay (90) 
= i=0) 
where 


Ber = Elz Pyy Phd (91) 


BuLion 0= 3537 — 4 


tJ 
t 


Ez; Px Ps = 0 (92) 


so that (90) may be reduced to 


[—] m 
21 Pri = 4AuP esi + whee a > loge (93) 
i=0 1=0 
where 
Aer = Elz Prep Pray il (94) 
and 
Aye = Elz, Pap Pir (95) 


The polynomials, P,.,,, may now be found from the lower order polynomuals 


i-—1 m 
a i 2 r 
Ar pay = 21 Px _ oe _ yy we tet (96) 
i=0 =O 
fata l 
P34) = (97) 
a WA gee al 


This completes the development of the bivariate Szego polynomials as pre- 
sented by Justice [Ref. 2]. The polynomial P,, will hold the coefficients for the k x/ 
prediction error filter. The desired causal prediction error filter, A,, 1s obtained by 
iotipiime Py bv 2z;*z;’ and dividing by p,{k,f to force the leading coefficient 
a,{0,0) = 1 as was assumed. The signal generator model is found from (30) with e(7,77,) 


taken as the input. 


l l 
—E—E ———E——eE (98) 


Alas aval 7! Bots! 
( 13 ) + Q\0-1 + Q91<2 +4321 <3 + aad 


This again assumes that A(z,, z,) 1s a causally invertible sequence. 


ha, 


D. LATTICE IMPLEMENTATION 
1. Analysis Lattice 

It is now a simple matter to utilize the orthogonal polynomials for lattice im- 
plementation. The first row is a function of z, only and the lattice structure is identical 
to the one-dimensional case. All subsequent polvnomials have a connection with each 
of the polynomials in the previous row and each of the polynomials in the same row but 
of lower order in z.. 

Using the familiar form of the Levinson algorithm the first row of the lattice 


filter may be drawn as shown in Figure 7, given by 
Agar = Aon — 22 44Aon (99) 
and 


Apiat = 22 Age — AaAon - (100) 


The orthogonal polynomials that have been developed will hold the filter coefficients in 
descending powers of z, andz,. That is, the lowest order term in P,,1s the term associated 
with the greatest delay in the filter difference equation. We could use (94), (95) and (96) 
to implement a first-quadrant filter that would be statistically equivalent to the third- 
quadrant causal filter we desire. If all the data are stored ahead of time, which 1s gen- 
erallv the case for 2-D sequences, this would not cause any computational difficulties. 
We will present the causal implementation of the lattice structure. The first quadrant 
development 1s similar. 

In order to correlate the powers of z,andz, with the delay and thus have a 
causal filter, we will define A,, = 2/*zz'P,, This will result in a realizable design and the 
form of the lattice filter. The first row of the filter is given bv (99) and (100). We now 
assume that row k has been constructed and that we proceed to rowkA+1. Multiplying 
(96) by Zz." 25* Vields 


i—1 m 

Peete eee! sen a) S(t), = Gees 

Ay) 22 Pays = 2 22 Pi — wie 22 Pegi — 21 22 DA ee Pee LOD) 
i=0 t=0 


that can be written in the form 





Figure 7. First row of the analysis lattice 
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Now with A4j,,, defined as 274257 Aerie 2) wesc 
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Equations (102) and (103) provide the means for implementing the two-dimensional 
lattice structure. Before proceeding to a schematic representation of the structure, it is 
worthwhile to explore the relationship of the J,,s and the filter coefficients. This will 
provide some insight into the operation of the lattice filter and make it apparent that the 
2's are sufficient to represent the filter. Assume that we have already found the lower 


order filters, Ago, Agi, Ag, and A, and now seek A,,. The coefficients that are available 


are shown in Figure 8. By (102) we find that 





Figure 8. The lower order polynomials for A,, 


=i s ay r ely r —ling r 
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oeemicer A,,(z,, z,) will then be of the form 


] 


£5) =) =e 
Ay1(2}5 22) = Aq +412, +4922 +4112) 2) (105) 


and (104) can now be used to equate the coefficients in A,, with the lower order 
polynomials. This is shown in Figure 9. The equations for the coefficients can be writ- 


eenedirectiy from Figure 9 as 


a, ,(00) = ao, (00) (106) 

a, (O01) = a,(01) — Agga@, (00) (107) 

a,,(10) = — £’o7a (92) (108) 

@ Ct) = — LAgoayo( 10) + 4’9,49,(01) + 2’o94@92(01) J (109) 
0 = ~L2' 9a 99(00) + 4’9141(00) + A%o242(00) J. (110) 


It can be seen that there are enough equations to determine all the filter coefficients of 
A,, Conversely, if the filters are known, then the /’s can be completely specified. This 
relationship allows the filter to be realized in either the direct form using the filter coef- 
Metents Or in the lattice structure using the /’s. If the same procedure is carried out 
Starting with (102), the coefficients in Ai, are found to correspond to those already found 
for A,,;. This is of course a requirement if the algorithm jis correct. 

The second row of the lattice structure described by (102) and (103) is shown in 
Figure 10 for a 2x3 filter. The first row of the filter is shown in Figure 7. All subsequent 
rows of the filter will have the same configuration. The structure can be seen to be quite 
similar to a one-dimensional lattice. There are. however, two important differences. 
First, there is a connection to all the lower-order predictors in the row. Also the inverse 
delay in the reverse error path appears to be a non-causal operator. These inverse delays 
fesleirom the definition of Aj, = z,"2;7A,Xz;", z') and force a specific ordering of points 
@em2 imputed to the lattice. For example, if we were computing the reverse error, 
é'(m1l,n2), for a filter with three horizontal filter coefficients, m= 2, then the z transforms 


@ieoi these errors are given by 


Eiji) = Ver 29)4}\(2) 22) (111) 


EA(lez2) ae Za) dance a aaa ar Bi us Era) (112) 
E\(z1,z2) = Y(z,, 2)(@goz1 23 me a2) ) (113) 


By referring to Figure 9 we see that the coefficients a,,(0,1) and a,,(0,0) must be 
in the same power of z,z, and likewise a,,(1,1) and a,(1,0) ._ Therefore, A;,(Z,, z,) as aoe 
multiphed by z, This would cause no problems if the direct form implementation was 


being used since 
2) Ejo(21,22) = Y(21,22)(Goo21 2, + M1022) - (114) 


There are no non-causal terms. In the lattice structure instead of linking actual filters 
we link error outputs. The calculation of ej,(m1,n2) requires ej,(71,n2 + 1). This requires 
the next horizontal value of e,, before the present value of e,, can be computed. The end 
result of this is that the previous row errors must be completely calculated and stored 
before the present row errors can be computed. This forces the input to proceed across 
rows. 

The filter shown in Figure 10 represents extending the filter order by rows. If 
at some point we Wish to extend the filter by columns the procedure 1s identical and the 
next column would have the same connections as the row extension with the subscripts 
transposed. 

2. Synthesis Lattice 

The all-pole model provides the means to synthesize the original signal from the 
error sequence. The synthesis filter can also be implemented in a lattice structure. 
Multiplving the recursion relationships (102) and (103) by Y(z) We arrive at the recursion 


for the error sequences. 
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Figure 9. Coefficient relationship between orthogonal polynomials. 
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Since Y(nl,n2) corresponds to e,, We must write (116) in decreasing order re- 
sulting in 
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The first row of this filter is again a function of one variable and is given by 
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Figure 10. One cell of a 2-D analysis lattice 


This 1s realized as shown in Figure 11, and the subsequent rows are given by 
(118). The signal flow graph is similar to Figure 10 with the signal flow inverted and the 
sign changes as indicated by (118). The node at e,, is shown in Figure 12. Computer 


simulation of these filters and the results obtained are presented in Chapter ITI. 


3 





Figure 11. First row of the synthesis lattice 
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Figure 12.  e, node of the synthesis lattice 
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Il. COMPUTER IMPLEMENTATION AND RESULTS 


A. IMPLEMENTATION 
1. Autocorrelation Function 
The autocorrelation function referred to throughout this thesis is actually an 
estimate based on a finite sample space. There are several ways to compute the esti- 
mated autocorrelation coefficients. For a large, complex-valued sample space the most 
efficient method is the use of Fourier transforms. For the small real-valued data space 
used in the examples of this study, the autocorrelation is estimated directly from the 


available data. The autocorrelation coefficients are given by 


Realn2) = > YK OL + hyn 4D) (121) 
k l 


2. Reflection Coefficients 
The basic algorithms for calculating the coefficients were presented by Justice 
[Ref. 2]. With some modification these algorithms are presented in this and the follow- 
ing section. 
The coefficients that link succesive order filters are derived in Chapter II. They 
will be referred to as reflection coefficients even though this term is descriptive of only 


the first row. The first row reflection coefficients were found to be 


Me >) Riera ; (122) 
k=0 


This result is simply the one dimensional Levinson algorithm appended with a second 
subscript. 
For subsequent rows there are two reflection coefficients, one for the forward 


polynomials and one for the reverse. Thev are obtained as 


Age = <2 Par Pray? (123) 


and 


a ae Be (124) 
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These are inner products over the autocorrelation space and are computed as 


l k+1 f 
kL >, ee es) (125) 
i. f=09s=) 
and 
m k l 
A kt = > Paik =jl> > Rae oll (126) 
J=0 r=0 s=0 


where any negative value in a subscript results in a zero. 
3. Filter Coefficients 
The first row of prediction error filters are polynomials of the form P,,. They 
can be calculated by the relationships 


ale 
Po.o(9,0) = (Ro) 2 


alls , 
Po,+191) = (1 = 14) *p3,(0,D) 2 (po (0 — 1) — Apo (9,)Po,0,/ — J) (128) 


sj =1 al 


(129) 
O</<m-] 


where mis the maximum filter order in the n2 direction and / 1s the filter order presently 
being computed. 
All subsequent rows will increase the filter order in the nl direction and can be 


recursively computed by 


Prat Md) = Pedi — Ly) - Saeed (ij) — S A Pe Ak — im — J). (130) 
I=} [=m—j 
Ui oaiael 
: (131) 
USGS 


This recursion can be carried on until the desired filter order is reached. The only 


polynomials required to increase the filter size are the preceding row and the lower order 


oo 


filters in the same row. It may be desired to normalize the filter equation, and this is 


easily carried out by 


Pres dl) 
Presi iy) = —a-—— (132) 
k+1,l 
where 
Ay) oa Press al (133) 
or 
k+1 1 
dy = [Pai hk + 1) > Rrercsta dr \TF | (134) 
f=0 s=0 
B. RESULTS 


J. System Model 
The algorithm which has been developed in this study is particularly useful in 
finding the autoregressive or all-pole model of the system which has generated these 
data. Several examples are presented in this section. Notice that the computer gener- 
ated results have the negative of the desired values. This results from our original deft- 


nition of the prediction of Yin (9) as 
Y == Lan Venln2 = Loan, Y ene eee (135) 
Figures 13 and 14 show various size system models that were derived from the 


data generated bv the difference equation 
YOul.2) = .4¥(n] — 1.n2) + .6Y¥(2l,n2—-—1) — .6¥(n] — ln2-—1) + X(n1,n2). (136) 


Figure 13 contains the system models resulting from driving the difference equation with 
an unit impulse. Figure 14 is the result of driving the system with a pseudo white noise 
sequence. The discrepancies in the filter coefficients can be mainly attributed to insuf- 
ficient ‘whiteness’ in the noise data. 
2. Spectral Estimation 
If we assume the signal, Y(n1,n2), is generated by the model of Figure 15, the 


power spectrum is given by 
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Figure 13. 
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System models for Eq.136 with an impulse input 


Figure 14. 
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System models for Eq.136 with a white noise sequence input 
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S(2), 22) = (FL Y(al,n2)])° (137) 


or 
S(z;, 29) = (FL B(l ,n2) * e(n1 n2) J) (138) 
or 
S(215 22) = (Ble1, 22) E(2), 22)" (139) 
If the input is a unit variance white noise sequence, the power spectrum be- 
comes 


S(24,22) = BY(21,29) (140) 


We can now use the all-pole system model to arrive at an estimate of the power spec- 


trum 


] | 


(5 | 14] 
(1,22) ( A(z, 29) A(z), 2,)A (2), 29) - 


where A(z,, z,) 1s the all-zero prediction error filter previously derived. 

We have considered three examples of two-dimensional all-pole power spectrum 
estimation. In all examples the data set size used to compute the autocorrelation func- 
tion is 32 x 32. In the first example a sequence consisting of two cosines at normalized 
frequencies (0.5,0.25) and (0.75,0.75) is corrupted by two-dimensional unit variance 
white noise. Figure 16 shows the surface plot of the estimated power spectrum with a 
4x 4 filter mask. The ideal response would be two impulse functions at the appropriate 
frequencies. The plot of Figure 16 shows several spurious peaks along a ridge which 1s 
illustrated by the contour plot in Figure 17. It can clearly be seen from the surface and 
contour plots in Figures 18 and 19 that the estimation of the power spectrum improves 
rapidly with an increased filter size. Figures 20 and 21 show the results using an § x § 
filter. The surface plot in Figure 20 clearly shows two discrete peaks in the power 
spectrum. Notice that there is a slight bias in the estimates. 

In the second example three cosines at frequencies (0.5,0.25), (0.75,0.75), and 
(0.25,0.75) are used to generate the data. The surface and contour plots of spectral es- 
timates are presented in Figures 22 through 27. The filter masks are square with sizes 


4x4,6x6,and8x&. The ability of these filters to distiguish the spectral peaks again 


oS 





Figure 15. Generating model for spectral estimates 


improves as the filter size increases from 4 x 4to 8 x 8. While the 4 x 4 filter, as shown 
in Figure 22, provides very little useful information, the 8 x 8 filter in Figure 26 clearly 
indicates three distinct peaks. 

The third example investigates the capability of the algorithm to differentiate 
between two closely spaced sinusoids at frequencies (0.4375,0.4375) and (0.5625,0.5625). 
As observed in Figures 28 through 31 the algorithm has not been able to distinguish 
between the two peaks. However, with a filter size 8 x 8 (see Figures 32 and 33), the 
spectral peaks become distinct even though there are some weak spurious peaks present. 

The examples which are presented in Figures 16 to 33 show the usefulness of the 
algorithm in estimation of the power spectrum. The ability of the algorithm to identify 
discrete peaks in the power density spectrum increases rapidly with filter order. As can 
be seen in Figures 32 and 33, even closely spaced peaks are readily identified by an 
eighth-order filter. This filter requires the computation of 64 coefficients which is much 
simpler than using Fourier transforms. All plots are scaled from 0 to | which represents 
the interval 0,z. 

3. Signal Synthesis 

The reflection coefficients for a single frequency data set were determined and 
used to implement a computer simulation of the synthesis lattice. The resulting lattice 
was then driven by a pseudo random white noise sequence. The spectruin of the output 
from the lattice is shown in Figure 34 for a 4 x 4 structure. Considering the small size 


of the filter and the poor quality of the input noise sequence, the results are acceptable. 
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Figure 16. Surface plot of spectral estimate: Y(nl,n2)=cos(nl 2/2 + n2 7/4) + 
Gesil 7 did + n2 7 3H) 
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Figure 17. Contour plot of spectral estimate: Y(nl,n2)=cos(ni 2/2 + n2 7/4) 
+ Cos(n lug 3/4 staeaeeeeo 4) 
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Figure 18. 
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Figure 19. Contour plot of spectral estimate: Y(ni,n2)=cos(nl 2/2 + n2 7/4) 
+ COs(h! a 3/4 4 Zen > 4) 
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Y(ni,n2)=cos(nl 2/2 + n2 7/4) + 


Surface plot of spectral estimate: 
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Figure 20. 
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Figure 21. Contour plot of spectral estimate: Y(nl,n2)=cos(nl 2/2 + n2 7/4) 
+ COs l x 3/4 “error 
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Figure 22. Surface plot of spectral estimate: Y(nil,n2)=cos(nl 2/2 + n2 2/4) + 
Coetal x 3/4 + n2 w 3,4) costal 2z/4yn2 x 3/4) 
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Figure 23. Contour plot of spectral estimate: Y(nl,n2)=cos(nl 2/2 + n2 7/4) 
+ cos(nl 7 3/4 -n2 3/4) cost ve oan 
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Figure 24. 
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Figure 25. Contour plot of spectral estimate: Y(nl,n2)=cos(nl 2/2 + n2 7/4) 


+ cos(nl 2 3/4 + n2 2 3/4) + cos(nl 2/4, n2 x 3/4) 
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Figure 26. 
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Figure 27. Contour plot of spectral estimate: Y(nl,n2)=cos(nl 2/2 + n2 7/4) 
+ cos(nl 2 3/4 + n2 7” 3/4) + cos(nl 2/4, n2 2 3/4) 
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Figure 28. Surface plot of spectral estimate: Y(nl,n2)=cos(nl 27/16 + n2 
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Figure 29, Contour plot of spectral estimate: Y(nI,n2)=cos(nl 2/7/16 + n2 
mi 16) cos(il a9) 16a 
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Figure 30. Surface plot of spectral estimate: Y(nl,n2)=cos(nl 27/16 + n2 
Moet cos(nlag 16+ n2 799716) 


38) 


6x6 FILTER 


92 AXIS 





fe) Or2 0.4 0.6 0.8 1.0 
61 AXIS 


Figure 31. Contour plot of spectral estimate: Y(ni,n2)=cos(nl 27/16 + n2 
RI/16) + cos(nlag IG - 2 a 
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Figure 32. Surface plot of spectral estimate: Y(nil,n2)=cos(nl 27/16 + n2 
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Figure 33. Contour plot of spectral estimate: Y(nl,n2)=cos(nl 27;16 + n2 
u7/16) + cos(niz9,16 + n2 2 9/16) 
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SYNTHESIZED SIGNAL 





Figure 34. Spectrum of signal generated from white noise: Y(n1,n2)=cos(nl 2/2 
seeilic 41} 


IV. CONCLUSIONS 


The linear prediction approach to analysis of stationary one-dimensional signals has 
direct applications in the study of two-dimensional signals. The solution of the 2-D 
normal equations in terms of orthogonal polynomials provides an efficient and effective 
method for system modeling and power spectral density estrmation. There also exist 
two-dimensional lattice structures that are analogous to their one-dimensional counter- 
parts. These lattices can be used as an alternate implementation of prediction error fil- 
ters and signal syntesis filters. 

The simple nature of the recursive solution to the 2-D normal equations developed 
in this thesis provides an excellent method for obtaining solutions to problems that are 
not dependent on the stability of the derived filter. For example, system modeling 
achieves very accurate results assuming the system being modeled 1s stable. Signal syn- 
thesis on the other hand can not be considered reliable since there is nothing to guar- 
antee a stable filter. This can be attributed to the lack of autocorrelation matching in 
two-dimensional signals. An excellent discussion of this property can be found 1n [Ref. 
4]. Simply stated, there are more independent autocorrelation coefficients than there are 
coefficients in the filters that we have found. We have ignored the second- and fourth- 
quadrant coeflicients which are independent from the first- and second-quadrant 
cocfliceints we have been using. 

Considering the computational expense of other spectrum estimation techniques the 
algorithm presented in this paper is valuable if used only for that purpose. The problem 
of signal synthesis, while more difficult, was met with some success as evidenced by 
Figure 34. This limited success invites further study into autoregressive modeling of 2-D 


sicnals. 
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